## January 2022.
## This file will create all figures and tables from "The Assessment Gap: Racial Inequalities in Property Taxation".

## Packages
library(data.table)
library(lfe)
library(stringr)
library(dplyr)
library(stargazer)
library(matlab)
library(ggplot2)

## Associated Function - Included
source("[PATH]/F - Get Number of Clusters for Stargazer.R")
source("[PATH]/F - Hedonic Helper Functions.R")
source("[PATH]/F - Distribution Helper Functions.R")
source("[PATH]/F - Helper Functions, Quantile Analysis.R")
times100=function(x) { 100*x }
options(scipen = 999)


#####################################
### Load Core Data ##################
#####################################

allgaps=fread("[PATH]/assessment gap, core dataset.csv",integer64 = "character") # Core dataset

allgaps$sr_property_id=as.character(allgaps$sr_property_id)  
allgaps=allgaps[setdiff(1:nrow(allgaps),which(allgaps$trans.prioryr==1)),] 
allgaps$logmkt=log(allgaps$sr_val_transfer)
allgaps$logassmt=log(allgaps$sa_val_assd)

gaps2=allgaps[allgaps$hmdasubset==1,c("sr_property_id","fips_10","group","year","geoid_10","state","tract","ag2","logmkt","logassmt","sr_date_transfer","sr_val_transfer",
                                      "black.factorvariable","blackhisp.factorvariable","allrace.factorvariable","other.factorvariable", "blackflag",  
                                      "black.acs","nonwhite.acs","blackhisp.acs","white.nh","pop.total.acs"
                                     )]
gaps2=gaps2[gaps2$year>=2005,]
gaps2=gaps2[complete.cases(gaps2),] 
gaps2$gtract=paste(gaps2$group,gaps2$tract,sep=".") 
gaps2$gtract=as.factor(gaps2$gtract)
gaps2$year.int=gaps2$year
gaps2$year=as.factor(gaps2$year)  
gaps2$group=as.factor(gaps2$group)  # Group is designation for a taxing jursidiction. (A "group" of overlapping governments.)

gaps2$black.factorvariable=as.factor(gaps2$black.factorvariable)
gaps2$blackhisp.factorvariable=as.factor(gaps2$blackhisp.factorvariable)
gaps2$allrace.factorvariable=as.factor(gaps2$allrace.factorvariable)
gaps2$other.factorvariable=as.factor(gaps2$other.factorvariable)

gaps2$black.factorvariable=relevel(gaps2$black.factorvariable,ref="W")
gaps2$blackhisp.factorvariable=relevel(gaps2$blackhisp.factorvariable,ref="W")
gaps2$allrace.factorvariable=relevel(gaps2$allrace.factorvariable,ref="W")
gaps2$other.factorvariable=relevel(gaps2$other.factorvariable,ref="W")


#####################################################
### HMDA Merge Analysis #############################
#####################################################

hmerge=allgaps[,c("sr_property_id","fips_10","group","year","geoid_10","state","tract","ag2","hmdasubset", "trans.prioryr",
                 "black.factorvariable","blackhisp.factorvariable","allrace.factorvariable","other.factorvariable", "blackflag",  
                 "black.acs","nonwhite.acs","blackhisp.acs","white.nh","pop.total.acs",
                 "ownerpercent","medianage.all","medianyear.owner","medianvalue","hhinc.med.total.acs",
                 "unemployment","notinlabforce","gini","hhfoodstamps"
                 ,"logmkt","logassmt" 
)]
hmerge=hmerge[hmerge$year>=2005,]     


hmerge$trans.prioryr[which(is.na(hmerge$trans.prioryr))]=0
hmerge$hmdasubset[which(is.na(hmerge$hmdasubset))]=0
hmerge$blackflag[which(is.na(hmerge$blackflag))]=9
hmerge$black.factorvariable[hmerge$black.factorvariable==""]="U"
hmerge$blackhisp.factorvariable[hmerge$blackhisp.factorvariable==""]="U"
hmerge$other.factorvariable[hmerge$other.factorvariable==""]="U"
hmerge=hmerge[complete.cases(hmerge),]


hmerge$gtract=paste(hmerge$group,hmerge$tract,sep=".") 
hmerge$gtract=as.factor(hmerge$gtract)
hmerge$year.int=hmerge$year 
hmerge$year=as.factor(hmerge$year)  
hmerge$group=as.factor(hmerge$group)

hmerge$black.factorvariable=as.factor(hmerge$black.factorvariable)
hmerge$blackhisp.factorvariable=as.factor(hmerge$blackhisp.factorvariable)
hmerge$allrace.factorvariable=as.factor(hmerge$allrace.factorvariable)
hmerge$other.factorvariable=as.factor(hmerge$other.factorvariable)

hmerge$black.factorvariable=relevel(hmerge$black.factorvariable,ref="W")
hmerge$blackhisp.factorvariable=relevel(hmerge$blackhisp.factorvariable,ref="W")
hmerge$allrace.factorvariable=relevel(hmerge$allrace.factorvariable,ref="W")
hmerge$other.factorvariable=relevel(hmerge$other.factorvariable,ref="W")

hmerge$medianvalue=log(as.numeric(hmerge$medianvalue))
hmerge$hhinc.med.total.acs=log(as.numeric(hmerge$hhinc.med.total.acs))
hmerge$pop.total.acs=log(hmerge$pop.total.acs)

######
### Bring in housing stock attributes
######

hedatt=fread("[PATH]/hedonic attributes, clean.csv",integer64 = "character")
hedatt$sa_property_id=as.character(hedatt$sa_property_id)

mergegaps=merge(hmerge,hedatt,by.x=c("sr_property_id","year.int"),by.y=c("sa_property_id","year"))  
mergegaps=mergegaps[setdiff(1:nrow(mergegaps),intersect(which(mergegaps$sa_nbr_bedrms==0),which(mergegaps$sa_nbr_bath==0))),] 
mergegaps=mergegaps[which(mergegaps$synthrooms>mergegaps$sa_nbr_stories),]
mergegaps$sa_fireplace_code=ifelse(mergegaps$sa_fireplace_code=="Y",1,0)
mergegaps$sa_pool_code=ifelse(mergegaps$sa_pool_code=="Y",1,0)
mergegaps$sa_patio_porch_code=ifelse(mergegaps$sa_patio_porch_code=="Y",1,0)
mergegaps=mergegaps[complete.cases(mergegaps),] 

for (i in 1:ncol(mergegaps)) {
  if (class(mergegaps[[i]])=="integer") { mergegaps[[i]]=as.numeric(mergegaps[[i]]) }
}

######
### Produce Estimates
######

testvars=c("black.acs","nonwhite.acs","blackhisp.acs","white.nh","pop.total.acs",
           "ownerpercent","medianage.all","medianyear.owner","medianvalue","hhinc.med.total.acs",
           "unemployment","notinlabforce","gini","hhfoodstamps",
           "logmkt","ag2",                                                                            
           "sa_sqft","bathadj", "sa_yr_blt","sa_patio_porch_code","sa_pool_code","sa_fireplace_code","sa_nbr_stories")

result=data.table(matrix(,nrow=length(testvars),ncol=2))
colnames(result)=c("coef.yeargroup","se.yeargroup")
result$variable=testvars

for (i in 1:length(testvars)) {
  var=testvars[i]
  zz3=summary(felm(formula(paste(var," ~ hmdasubset | group:year | 0 | group",sep="")), data=mergegaps))
  
  result$coef.yeargroup[result$variable==var]=zz3$coefficients[1,1]
  result$se.yeargroup[result$variable==var]=zz3$coefficients[1,2]
}

## Order based on order of test vars (and hence result output) - DO NOT CHANGE. 
result$label=c("Population Share Black","Population Share Non-White","Population Share Black or Hispanic","Population Share White","Population (log)",
               "Owner Percentage","Median Age (Yrs)","Median Year Purchased","Median Home Value (log)","Median HH Income (log)","Unemployment Rate",
               "Not In Labor Force Share","Gini Coefficient","Share SNAP","Transaction Price (log)","Assessment Ratio (log)","Square Feet","Nbr Baths","Year Built","Patio or Porch (Binary)",
               "Pool (Binary)","Fireplace (Binary)","Number of Stories")

table=rbind(result[,c("label","coef.yeargroup")],result[,c("label","se.yeargroup")],fill=T)
table=table[order(match(table$label,result$label)),]
table$output=ifelse(!is.na(table$coef.yeargroup),table$coef.yeargroup,table$se.yeargroup)
table$output=round(table$output,5)
table$label[seq(2,nrow(table),2)]=""
table$output[seq(2,nrow(table),2)]=paste("(",table$output[seq(2,nrow(table),2)],")",sep="")
colnames(table)[which(colnames(table)=="label")]="Tract or Property Attribute"
colnames(table)[which(colnames(table)=="output")]="Difference"

## Add stars
table$tstat=NA
table$tstat[seq(1,nrow(table),2)]=table$coef.yeargroup[seq(1,nrow(table),2)]/table$se.yeargroup[seq(2,nrow(table),2)]
table$stars=ifelse(abs(table$tstat)>2.5758,"***",ifelse(abs(table$tstat)>1.959964,"**",ifelse(abs(table$tstat)>1.64485,"*","")))
table$Difference[seq(1,nrow(table),2)]=paste(table$Difference,table$stars,sep="")[seq(1,nrow(table),2)]


## Output Table
stargazer(table[,c("Tract or Property Attribute","Difference")],type="latex",
          summary=F,align=T,float=F,rownames = F,
          out = NULL #"[OUTPUT PATH]"
)

## Output Coefficients
fwrite(result,"[OUTPUT PATH]")

#####################################################
### Baseline AG Result ##############################
#####################################################

## "ag2" is logged assessment ratio: ln(A/M)

race.jur1=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=gaps2)
race.jur2=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group,data=gaps2)
race.jur3=felm(ag2 ~ other.factorvariable | group:year | 0 | group, data=gaps2)

########
### At tract level 
#######

race.tract1=felm(ag2 ~  black.factorvariable | gtract:year | 0 | group,data=gaps2)
race.tract2=felm(ag2 ~ blackhisp.factorvariable | gtract:year | 0 | group,data=gaps2)
race.tract3=felm(ag2 ~ other.factorvariable | gtract:year | 0 | group, data=gaps2)

########
### At block-group level
########

## Load block group identifiers
bgs=fread("[PATH]/block group assignments.csv",colClasses = "character",
          select = c("sa_property_id","BG.ID"))

gaps2=merge(gaps2,bgs,by.x="sr_property_id",by.y="sa_property_id",all.x=T,all.y=F)
gaps2$BG.ID[which(nchar(gaps2$BG.ID)==0)]=NA
gaps2$BG.ID=str_pad(gaps2$BG.ID,width = 12,pad = "0",side = "left")

gaps2$govbg=paste(gaps2$group,gaps2$BG.ID,sep=".")
gaps2$govbg=factor(gaps2$govbg)

race.bg1=felm(ag2 ~  black.factorvariable | govbg:year | 0 | group,data=gaps2)
race.bg2=felm(ag2 ~ blackhisp.factorvariable | govbg:year | 0 | group,data=gaps2)
race.bg3=felm(ag2 ~ other.factorvariable | govbg:year | 0 | group,data=gaps2)


########
### Produce Table
########

baselineB=list(race.jur1,race.tract1,race.bg1)
baselineBclust=getNclust(baselineB)
stargazer(baselineB,type="latex",
          title="Assessment Gap, Black Homeowners",
          keep=c("black.factorvariableB"),
          covariate.labels = c("Black Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects","Jurisd-Year","Tract-Year","BG-Year"),c("No. Clusters",baselineBclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

baselineBH=list(race.jur2,race.tract2,race.bg2)
baselineBHclust=getNclust(baselineBH)
stargazer(baselineBH,type="latex",
          title="Assessment Gap, Black or Hispanic Homeowners",
          keep=c("blackhisp.factorvariableBH"),
          covariate.labels = c("Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(baselineB))),c("No. Clusters",baselineBclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

### Appendix Table 
other=list(race.jur3,race.tract3,race.bg3)
otherclust=getNclust(other)
stargazer(other,type="latex",
          title="Inequality for non-Black, non-Hispanic Minorities",
          keep=c("other.factorvariableO"),
          covariate.labels = c("Other Nonwhite Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects","Jurisd-Yr","Tract-Yr","BG-Yr"),c("No. Clusters",otherclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

#####################################################
### State Distribution ##############################
#####################################################

########
### Black Homeowners
########

state.base=distrib(gaps2,"state",level = "group")
state.base=state.base[complete.cases(state.base),]
state.base=state.base[order(state.base$coef,decreasing = T),]
state.base$xaxis=1:nrow(state.base)
rows=nrow(state.base)

pdf("[OUTPUT PATH]",paper="a4r",width=10.5,height = 8)
plot(state.base$xaxis,state.base$coef,ylim=c(-.15,.45),type="p",pch=16, main="Assessment Gap, by State",ylab="Assmt Gap, Black Residents",xaxt="n",xlab="",cex=1.6,cex.lab=1.5,cex.main=1.5)
par(new=T)
arrows(1:rows,state.base$lb,1:rows,state.base$ub,angle=90,code=3,length=.02,lwd=2)
abline(h=0)
endpoint=rows*1.0
text(x=seq(.65,endpoint,by=(endpoint-.65)/(rows-1)),y=-.12,pos=3,par("usr")[3]-.25,srt=60, adj=1,xpd=T,labels=state.base$state,cex=.95)
dev.off()

########
### Black or Hispanic Homeowners
########

state.base.bh=distrib.bh(gaps2,"state",level = "group")
state.base.bh=state.base.bh[complete.cases(state.base.bh),]
state.base.bh=state.base.bh[order(state.base.bh$coef,decreasing = T),]
state.base.bh$xaxis=1:nrow(state.base.bh)
rows=nrow(state.base.bh)

pdf("[OUTPUT PATH]",paper="a4r",width=10.5,height = 8)
plot(state.base.bh$xaxis,state.base.bh$coef,ylim=c(-.15,.45),type="p",pch=16, main="Assessment Gap, by State",ylab="Assmt Gap, Black or Hispanic Residents",xaxt="n",xlab="",cex=1.6,cex.main=1.5,cex.lab=1.5)
par(new=T)
arrows(1:rows,state.base.bh$lb,1:rows,state.base.bh$ub,angle=90,code=3,length=.02,lwd=2)
abline(h=0)
endpoint=rows*1.0
text(x=seq(.65,endpoint,by=(endpoint-.65)/(rows-1)),y=-.12,pos=3,par("usr")[3]-.25,srt=60, adj=1,xpd=T,labels=state.base.bh$state,cex=.95)
dev.off()


#####################################################
### County Distribution #############################
#####################################################

countydistrib=distrib.county(gaps2,"fips_10") 
countydistrib=countydistrib[which(!is.na(countydistrib$coef)),]
countydistrib=countydistrib[order(countydistrib$coef,decreasing = T),]
countydistrib$xaxis=1:nrow(countydistrib)

pdf("[OUTPUT PATH]",paper="a4r",width=10.5,height = 8)
plot(countydistrib$xaxis,countydistrib$coef,ylim=c(-.2,.5),type="p",pch=16,cex=1.6,main="County Level Distribution",xaxt="no",xlab=paste(nrow(countydistrib)," Counties Total"),ylab="Assmt Gap, Black Resident",cex.lab=1.5)
abline(h=0)
dev.off()

#####################################################
### Conduct Regressions by County ###################
#####################################################

gaps2$fips=substr(gaps2$tract,1,5)
gaps2$fips=as.factor(gaps2$fips)

county1=felm(ag2 ~  black.factorvariable | fips:year | 0 | fips,data=gaps2)
county2=felm(ag2 ~ blackhisp.factorvariable | fips:year | 0 | fips,data=gaps2)
county3=felm(ag2 ~ other.factorvariable | fips:year | 0 | fips, data=gaps2)

county=list(county1,county2,county3)
countyclust=getNclust(county)

stargazer(county,type="latex",
          title="Assessment Gap Within Counties Instead of Jurisdictions",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder", "Other Nonwhite Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("County-Year",3)),c("No. Clusters",countyclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

#####################################################
### By Year #########################################
#####################################################

B2005=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2005,])
B2006=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2006,])
B2007=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2007,])
B2008=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2008,])
B2009=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2009,])
B2010=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2010,])
B2011=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2011,])
B2012=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2012,])
B2013=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2013,])
B2014=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2014,])
B2015=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2015,])
B2016=felm(ag2 ~ black.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2016,])

BH2005=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2005,])
BH2006=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2006,])
BH2007=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2007,])
BH2008=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2008,])
BH2009=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2009,])
BH2010=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2010,])
BH2011=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2011,])
BH2012=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2012,])
BH2013=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2013,])
BH2014=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2014,])
BH2015=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2015,])
BH2016=felm(ag2 ~ blackhisp.factorvariable | group | 0 | group,data=gaps2[gaps2$year.int==2016,])

Byear=list(B2005,B2006,B2007,B2008,B2009,B2010,B2011,B2012,B2013,B2014,B2015,B2016)
Byearclust=getNclust(Byear)
stargazer(Byear,type="latex",
          title="Assessment Gap by Year",
          keep=c("black.factorvariableB"),
          covariate.labels = c("Black Mortgage Holder"),
          column.labels = as.character(c(2005:2016)),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(Byear))),c("No. Clusters",Byearclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

BHyear=list(BH2005,BH2006,BH2007,BH2008,BH2009,BH2010,BH2011,BH2012,BH2013,BH2014,BH2015,BH2016)
BHyearclust=getNclust(BHyear)
stargazer(BHyear,type="latex",
          title="Assessment Gap by Year",
          keep=c("blackhisp.factorvariableBH"),
          covariate.labels = c("Black or Hispanic Mortgage Holder"),
          column.labels = as.character(2005:2016),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(BHyear))),c("No. Clusters",BHyearclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

#############################################
### By Segregation Index ####################
#############################################

seggaps=gaps2

disind=fread("[PATH]/index of dissimilarity.csv") ## Index of Dissimilarity created from 2000 Census data (tract level).
disind=disind[intersect(which(!is.na(disind$wb)),which(!is.na(disind$wbh))),] 
disind[,"quintileB" := cut(disind$wb,quantile(disind$wb,seq(0,1,1/5)),include.lowest=T,labels=LETTERS[1:5])]
disind[,"quintileBH" := cut(disind$wbh,quantile(disind$wb,seq(0,1,1/5)),include.lowest=T,labels=LETTERS[1:5])]
disind[,"decileB" := cut(disind$wb,quantile(disind$wb,seq(0,1,1/10)),include.lowest=T,labels=LETTERS[1:10])]
disind[,"decileBH" := cut(disind$wbh,quantile(disind$wbh,seq(0,1,1/10)),include.lowest=T,labels=LETTERS[1:10])]

oldfips=allgaps[allgaps$hmdasubset==1,c("sr_property_id","fips_00")]
oldfips=oldfips[!duplicated(oldfips),]

seggaps=merge(seggaps,oldfips,by.x="sr_property_id",by.y="sr_property_id",all.x=T,all.y=F)
seggaps=merge(seggaps,disind,by.x="fips_00",by.y="fips")

seg1=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$quintileB=="A",])
seg2=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$quintileB=="B",])
seg3=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$quintileB=="C",])
seg4=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$quintileB=="D",])
seg5=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$quintileB=="E",])


seg6=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="A",])
seg7=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="B",])
seg8=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="C",])
seg9=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="D",])
seg10=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="E",])
seg11=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="F",])
seg12=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="G",])
seg13=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="H",])
seg14=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="I",])
seg15=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileB=="J",])


seg16=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="A",])
seg17=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="B",])
seg18=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="C",])
seg19=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="D",])
seg20=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="E",])
seg21=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="F",])
seg22=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="G",])
seg23=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="H",])
seg24=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="I",])
seg25=felm(ag2 ~  blackhisp.factorvariable | group:year | 0 | group,data=seggaps[seggaps$decileBH=="J",])



segdecB=list(seg6,seg7,seg8,seg9,seg10,seg11,seg12,seg13,seg14,seg15)
segdecBclust=getNclust(segdecB)
stargazer(segdecB,type="latex",
          keep=c("black.factorvariableB"),
          covariate.labels = c("Black Mortgage Holder"),
           keep.stat = c("n","rsq"),
           dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jur-Yr",length(segdecB))),c("No. Clusters",segdecBclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

segdecBH=list(seg16,seg17,seg18,seg19,seg20,seg21,seg22,seg23,seg24,seg25)
segdecBHclust=getNclust(segdecBH)
stargazer(segdecBH,type="latex",
          keep=c("blackhisp.factorvariableBH"),
          covariate.labels = c("Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jur-Yr",length(segdecB))),c("No. Clusters",segdecBHclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

segresults=data.table(matrix(0,ncol=5,nrow=10))
colnames(segresults)=c("decile","coef.B","se.B","coef.BH","se.BH")
for (i in 1:10) {
  segresults$decile[i]=i
  bb=segdecB[[i]]
  segresults$coef.B[i]=bb$coefficients[which(rownames(bb$coefficients)=="black.factorvariableB")]
  segresults$se.B[i]=bb$cse[which(names(bb$cse)=="black.factorvariableB")]
  hh=segdecBH[[i]]
  segresults$coef.BH[i]=hh$coefficients[which(rownames(hh$coefficients)=="blackhisp.factorvariableBH")]
  segresults$se.BH[i]=hh$cse[which(names(hh$cse)=="blackhisp.factorvariableBH")]
}

segresults$decile=as.factor(segresults$decile)

## Output Figure
segB=segresults[,c("coef.B","se.B")]
colnames(segB)=c("coef","se")
segB$lb=segB$coef-1.96*segB$se
segB$ub=segB$coef+1.96*segB$se
pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(segB,names=c("Low",2:9,"High"),clr="blue",ylim=c(0,.25),num=10,ylab="Assessment Gap, Black Homeowners",main="Assessment Gap by Segregation Decile")
dev.off()

segBH=segresults[,c("coef.BH","se.BH")]
colnames(segBH)=c("coef","se")
segBH$lb=segBH$coef-1.96*segBH$se
segBH$ub=segBH$coef+1.96*segBH$se
pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(segBH,names=c("Low",2:9,"High"),clr="blue",ylim=c(0,.25),num=10,ylab="Assessment Gap, Black or Hispanic Homeowners",main="Assessment Gap by Segregation Decile")
dev.off()


#############################################
### Carveout for Cook County ################
#############################################

cookgaps=gaps2[gaps2$fips_10 == 17031,] 

cookB=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=cookgaps)
cookB.tract=felm(ag2 ~  black.factorvariable | gtract:year | 0 | group,data=cookgaps)
cookB.bg=felm(ag2 ~  black.factorvariable | govbg:year | 0 | group,data=cookgaps)

cookBH=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group,data=cookgaps)
cookBH.tract=felm(ag2 ~ blackhisp.factorvariable | gtract:year | 0 | group,data=cookgaps)
cookBH.bg=felm(ag2 ~ blackhisp.factorvariable | govbg:year | 0 | group,data=cookgaps)

cook=list(cookB,cookB.tract,cookB.bg,cookBH,cookBH.tract,cookBH.bg)
cookclust=getNclust(cook)

stargazer(cook,type="latex",
          title="Inequality in Cook County",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects","Jur-Yr","Tract-Yr","BG-Yr","Jur-Yr","Tract-Yr","BG-Yr"),c("No. Clusters",cookclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

#######################################################################
### Neighborhood Composition Regs with HMDA and ACS, no Controls ######
#######################################################################

gaps.nc=allgaps[allgaps$hmdasubset==1,c("black.acs","non.wbh","fips_10","group","blackhisp.acs","year","geoid_10","state","tract","ag2",
                                        "hhinc.med.total.acs", "ownerpercent","medianage.all","medianyear.owner","medianvalue",  
                                        "unemployment","notinlabforce","gini","hhfoodstamps",  
                                        "black.factorvariable","blackhisp.factorvariable","other.factorvariable"
)]

gaps.nc$hhinc.med.total.acs=as.numeric(gaps.nc$hhinc.med.total.acs) 
gaps.nc$medianvalue=as.numeric(gaps.nc$medianvalue)
gaps.nc$hhinc.med.total.acs=log(gaps.nc$hhinc.med.total.acs)
gaps.nc$medianvalue=log(gaps.nc$medianvalue)


gaps.nc=gaps.nc[complete.cases(gaps.nc),]
gaps.nc$year=as.factor(gaps.nc$year)
gaps.nc$group=as.factor(gaps.nc$group)

gaps.nc$black.factorvariable=as.factor(gaps.nc$black.factorvariable)
gaps.nc$blackhisp.factorvariable=as.factor(gaps.nc$blackhisp.factorvariable)
gaps.nc$other.factorvariable=as.factor(gaps.nc$other.factorvariable)

gaps.nc$black.factorvariable=relevel(gaps.nc$black.factorvariable,ref="W")
gaps.nc$blackhisp.factorvariable=relevel(gaps.nc$blackhisp.factorvariable,ref="W")
gaps.nc$other.factorvariable=relevel(gaps.nc$other.factorvariable,ref="W")

demo1=felm(ag2 ~ black.acs + black.factorvariable | group:year | 0 | group,data=gaps.nc)
demo2=felm(ag2 ~ blackhisp.acs + blackhisp.factorvariable | group:year | 0 | group,data=gaps.nc)

demoregs.main=list(demo1, demo2)
demoregsclust.main=getNclust(demoregs.main)

stargazer(demoregs.main,type="latex",
          title="",
          keep=c("black.acs","black.factorvariableB","blackhisp.acs","blackhisp.factorvariableBH"),
          order=c("black.factorvariableB","black.acs","blackhisp.factorvariableBH","blackhisp.acs"),
          covariate.labels = c("Black Mortgage Holder","Black Share","Black or Hispanic Mortgage Holder","Black or Hispanic Share"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(demoregs.main))),c("No. Clusters",demoregsclust.main)),
          digits=3,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          column.sep.width="10pt",
          out = NULL #"[OUTPUT PATH]"
)


########
### All controls - in standard deviations
########

for (i in c("black.acs","blackhisp.acs","non.wbh","hhinc.med.total.acs", "ownerpercent", "medianvalue", "medianage.all", 
            "unemployment", "notinlabforce", "gini", "hhfoodstamps")) {
  newcol=paste(i,".sd",sep="")
  gaps.nc[[newcol]]=gaps.nc[[i]]/sd(gaps.nc[[i]])
  
}

both1b=felm(ag2 ~ black.acs.sd + black.factorvariable + 
              hhinc.med.total.acs.sd + ownerpercent.sd + 
              unemployment.sd + gini.sd + hhfoodstamps.sd + medianage.all.sd | group:year | 0 | group,data=gaps.nc)

both2b=felm(ag2 ~ blackhisp.acs.sd + blackhisp.factorvariable + 
              hhinc.med.total.acs.sd + ownerpercent.sd + 
              unemployment.sd + gini.sd + hhfoodstamps.sd + medianage.all.sd | group:year | 0 | group,data=gaps.nc)

gapsregs2b=list(both1b, both2b)
clusterlength2b=getNclust(gapsregs2b)
stargazer(gapsregs2b,type="latex",
          title="gaps with controls",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","black.acs","blackhisp.acs",
                 "hhinc.med.total.acs","ownerpercent","unemployment", "hhfoodstamps", "gini","medianage.all"),
          order=c("black.factorvariable","black.acs","blackhisp.factorvariable","blackhisp.acs",
                  "hhinc.med.total.acs","unemployment","hhfoodstamps","ownerpercent","gini","medianage.all"),
          covariate.labels = c("Black Mortgage Holder","Black Share","Black or Hispanic Mortgage Holder","Black or Hispanic Share",
                               "Median HH Income","Unemployment","SNAP Assistance","Owner Percentage", "GINI Coef","Median Age"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(gapsregs2b))),
                         c("No. Clusters",clusterlength2b)),
          digits=3,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          column.sep.width="10pt",
          out = NULL #"[OUTPUT PATH]"
)

###############################################
### Sample Split on Cap Exists/Binds ##########
###############################################

capgaps=gaps2
capgaps=capgaps[capgaps$year.int>=2006,] # Lincoln Inst data from 2006.

######
### Load Data
######

caps=fread("[PATH]/capsdata.csv")

hpis=fread("[PATH]/clean annual zip code hpis (fhfa and zillow).csv")
hpis$zipcode=str_pad(hpis$zipcode,5,pad="0",side="left")

allzips=fread("[PATH]/zip codes for clean AG sample.csv",colClasses = "chr")
allzips$arc_zip=str_pad(allzips$arc_zip,5,pad="0",side="left")

######
### Analysis
######

statecaps=caps[caps$state %in% unique(caps$state[caps$statewide_policy==1 & caps$capdummy==1]),] 

capgaps=merge(capgaps,statecaps[,c("state","year","capdummy")],by.x=c("state","year.int"),by.y=c("state","year"),all.x=T,all.y=F) 
capgaps$capdummy[which(is.na(capgaps$capdummy))]=0 

capgaps=merge(capgaps,statecaps[,c("state","year","growthcutrate")],by.x=c("state","year.int"),by.y=c("state","year"),all.x=T,all.y=F) 

oneoff.cook=17031
oneoff.ny=setdiff(as.numeric(trimws(unlist(strsplit(paste(setdiff(unique(caps$fipscode),""),collapse=", "),",")))),oneoff.cook)
capgaps$capdummy[capgaps$fips_10 %in% oneoff.cook & capgaps$year.int<2015]=0  
capgaps$capdummy[capgaps$fips_10 %in% oneoff.ny]=1 
capgaps$growthcutrate[capgaps$fips_10 %in% oneoff.ny]=.06

capgaps=merge(capgaps,allzips,by.x="sr_property_id",by.y="sa_property_id",all.x=T,all.y=F) 
capgaps=merge(capgaps,hpis[,c("zipcode","year","hpi.fhfapref","hpi.zillpref")],by.x=c("arc_zip","year.int"),by.y=c("zipcode","year"),all.x=T,all.y=F)

capgaps$capbinds1=ifelse(capgaps$hpi.fhfapref<capgaps$growthcutrate,0,1)
capgaps$capbinds2=ifelse(capgaps$hpi.zillpref<capgaps$growthcutrate,0,1)


#######
### Regressions
#######

cap1=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==0])
cap2=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1])
#cap3=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1 & capgaps$capbinds1==1])  # No meaningful difference
cap3a=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1 & capgaps$capbinds2==1])

cap4=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==0])
cap5=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1])
#cap6=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1 & capgaps$capbinds1==1])  # No meaningful difference.
cap6a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=capgaps[capgaps$capdummy==1 & capgaps$capbinds2==1])

cap=list(cap1,cap2,cap3a,cap4,cap5,cap6a) 
capclust=getNclust(cap)
stargazer(cap,type="latex",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          column.labels = c("No Cap","Cap Exists","Cap Exists and Binds","No Cap","Cap Exists","Cap Exists and Binds"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(cap))),c("No. Clusters",capclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

###############################################
### Sample Split on Re-evaluation Freq ########
###############################################

cycles=fread("[PATH]/reeval cycles.csv")
colnames(cycles)[which(colnames(cycles)=="Physical")]="physicaleval"

cyclegaps=gaps2
cyclegaps=cyclegaps[cyclegaps$year.int>=2006,]  # Lincoln Institute Data starts 2006
cyclegaps=merge(cyclegaps,cycles[,c("State","Year","Cycle","physicaleval")],by.x=c("state","year.int"),by.y=c("State","Year"),all.x=T,all.y=F)
cyclegaps$Cycle[which(cyclegaps$fips_10==17031)]=3 

cycle1=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==1])
cycle2=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==2])
cycle3=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==3])
cycle4=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==4])
cycle5=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==5])
cycle6=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==6])
cycle8=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==8])
cycle9=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==9])
cyclenone=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==99])

cycle=list(cycle1,cycle2,cycle3,cycle4,cycle5,cycle6,cycle8,cycle9,cyclenone)
cycleclust=getNclust(cycle)
stargazer(cycle,type="latex",
          keep=c("black.factorvariableB"),
          covariate.labels = c("Black Mortgage Holder"),
          column.labels = c("1yr","2yrs","3yrs","4yrs","5yrs","6yrs","8yrs","9yrs","none"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(cycle))),c("No. Clusters",cycleclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

cycle1a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==1])
cycle2a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==2])
cycle3a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==3])
cycle4a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==4])
cycle5a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==5])
cycle6a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==6])
cycle8a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==8])
cycle9a=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==9])
cyclenonea=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=cyclegaps[cyclegaps$Cycle==99])

cyclea=list(cycle1a,cycle2a,cycle3a,cycle4a,cycle5a,cycle6a,cycle8a,cycle9a,cyclenonea)
cycleclusta=getNclust(cyclea)
stargazer(cyclea,type="latex",
          keep=c("blackhisp.factorvariableBH"),
          covariate.labels = c("Black or Hispanic Mortgage Holder"),
          column.labels = c("1yr","2yrs","3yrs","4yrs","5yrs","6yrs","8yrs","9yrs","none"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",length(cyclea))),c("No. Clusters",cycleclusta)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

############################################
### Load and Prepare Hedonic File ##########
############################################

hedgaps=merge(allgaps[allgaps$hmdasubset==1,],hedatt,by.x=c("sr_property_id","year"),by.y=c("sa_property_id","year"))  
hedgaps=hedgaps[hedgaps$year>=2005,] 
hedgaps=hedgaps[setdiff(1:nrow(hedgaps),intersect(which(hedgaps$sa_nbr_bedrms==0),which(hedgaps$sa_nbr_bath==0))),] 
hedgaps=hedgaps[which(hedgaps$synthrooms>hedgaps$sa_nbr_stories),] 

racevars=c("blackhisp.acs","black.acs","non.wbh")
neighborhoodvars=c("hhinc.med.total.acs","ownerpercent","unemployment", "hhfoodstamps", "gini","medianage.all")
hedvars=c("sa_sqft","betterroom","betterbath","sa_nbr_stories","sa_patio_porch_code","sa_pool_code","sa_fireplace_code","sa_yr_blt","synthrooms","bathadj")
usevars=c("sa_val_assd", "sr_val_transfer","tract", "year","group",racevars,neighborhoodvars,hedvars)
hedregfile=hedgaps[,usevars,with=F]

hedregfile$betterbath=as.factor(hedregfile$betterbath)
hedregfile=hedregfile[setdiff(1:nrow(hedregfile),which(hedregfile$betterbath=="")),]
hedregfile$betterbath=relevel(hedregfile$betterbath,ref="[0,1]")
hedregfile$betterroom=as.factor(hedregfile$betterroom)
hedregfile$betterroom=relevel(hedregfile$betterroom,ref="[0,1]")
hedregfile$sa_nbr_stories=as.factor(hedregfile$sa_nbr_stories)

hedregfile$hhinc.med.total.acs=log(as.numeric(hedregfile$hhinc.med.total.acs))
hedregfile$logassmt=log(hedregfile$sa_val_assd)
hedregfile$logmkt=log(hedregfile$sr_val_transfer)

hedregfile=hedregfile[complete.cases(hedregfile),]
hedregfile$tract=as.factor(hedregfile$tract)
hedregfile$group=as.factor(hedregfile$group)
hedregfile$year=as.factor(hedregfile$year)

for (i in 1:ncol(hedregfile)) {
  if (class(hedregfile[[i]])=="integer") { hedregfile[[i]]=as.numeric(hedregfile[[i]]) }
}

####################################################
### Hedonic Regression Table and Bar Chart #########
####################################################

hedregfile$sa_patio_porch_code=ifelse(hedregfile$sa_patio_porch_code=="Y",1,0)
hedregfile$sa_fireplace_code=ifelse(hedregfile$sa_fireplace_code=="Y",1,0)
hedregfile$sa_pool_code=ifelse(hedregfile$sa_pool_code=="Y",1,0)

for (i in c(racevars,neighborhoodvars,c("sa_sqft","bathadj", "sa_yr_blt","sa_patio_porch_code","sa_pool_code","sa_fireplace_code"))) {
  newcol=paste(i,"sd",sep=".")
  hedregfile[[newcol]]=hedregfile[[i]]/sd(hedregfile[[i]])
}

## NOTE: The helper function used to produce barplot data relies on variable order written below. DON'T CHANGE!
a1sd=felm(logassmt ~ black.acs.sd +
            hhinc.med.total.acs.sd + ownerpercent.sd + unemployment.sd + gini.sd + hhfoodstamps.sd +
            sa_sqft.sd +  bathadj.sd + sa_patio_porch_code.sd + sa_pool_code.sd + sa_fireplace_code.sd + sa_yr_blt.sd |
            group:year | 0 | group, data=hedregfile) 

m1sd=felm(logmkt ~ black.acs.sd +
            hhinc.med.total.acs.sd + ownerpercent.sd + unemployment.sd + gini.sd + hhfoodstamps.sd +
            sa_sqft.sd +  bathadj.sd + sa_patio_porch_code.sd + sa_pool_code.sd + sa_fireplace_code.sd + sa_yr_blt.sd |
            group:year | 0 | group, data=hedregfile) 


cont.sd=compare(a1sd,m1sd)
cont.sd
pdf("[OUTPUT PATH]",paper="a4r",width=10.5,height = 8)
print(getbarplot(cont.sd,"Implied Elasticity of Assessment Ratio to 1 SD Shift","Black Share",vers = "diff"))
dev.off()

## NOTE: The helper function used to produce barplot data relies on variable order written below. DON'T CHANGE!
a2sd=felm(logassmt ~ blackhisp.acs.sd +
            hhinc.med.total.acs.sd + ownerpercent.sd + unemployment.sd + gini.sd + hhfoodstamps.sd +
            sa_sqft.sd +  bathadj.sd + sa_patio_porch_code.sd + sa_pool_code.sd + sa_fireplace_code.sd + sa_yr_blt.sd |
            group:year | 0 | group, data=hedregfile) 

m2sd=felm(logmkt ~  blackhisp.acs.sd +
            hhinc.med.total.acs.sd + ownerpercent.sd + unemployment.sd + gini.sd + hhfoodstamps.sd + 
            sa_sqft.sd +  bathadj.sd + sa_patio_porch_code.sd + sa_pool_code.sd + sa_fireplace_code.sd + sa_yr_blt.sd |
            group:year | 0 | group, data=hedregfile) 

contbh.sd=compare(a2sd,m2sd)
pdf("[OUTPUT PATH]",paper="a4r",width=10.5,height = 8)
print(getbarplot(contbh.sd,"Implied Elasticity of Assessment Ratio to 1 SD Shift","Black/Hispanic Share",vers="diff",spread=.44))
dev.off()

pdf("[OUTPUT PATH]",paper="a4r",width=8,height = 10.5)
print(getbarplot(contbh.sd,"Implied Elasticity of Assessment Ratio to 1 SD Shift","Black/Hispanic Share",vers="diff",spread=.44))
dev.off()


########
##### Produce Table
########

allhed.main=list(m1sd,a1sd,m2sd,a2sd)
allhedclust.main=getNclust(allhed.main)

stargazer(allhed.main,type="latex",
          title="Hedonic Regressions",
          keep=c("black.acs.sd","blackhisp.acs",
                 "hhinc","unemployment","hhfoodstamps","ownerpercent.sd","gini.sd",
                 "sa_sqft","bathadj.sd","sa_yr_blt.sd"),
          order=c("black.acs.sd","blackhisp.acs",
                  "hhinc","unemployment","hhfoodstamps","ownerpercent.sd","gini.sd",
                  "sa_sqft","bathadj.sd","sa_yr_blt.sd"),
          covariate.labels = c("Black Share","Black or Hispanic Share",
                               "Median HH Income","Unemployment","SNAP Share","Owner Share","GINI",
                               "Square Feet","Bathrooms","Year Built"),
          keep.stat = c("n","rsq"),
          dep.var.labels = c("Market","Assessment","Market","Assessment"),
          add.lines=list(c("Other Attributes",rep("Y",length(allhed.main))),c("Fixed Effects",rep("Jurisd-Year",length(allhed.main))),c("No. Clusters",allhedclust.main)),
          digits=3,digits.extra = 1,float=F,dep.var.caption = "",
          notes = c("Not shown: coefficients on indicators for","patio, pool, and fireplace"),
          out = NULL #"[OUTPUT PATH]"
)

################################################
### Price Bins and Synthetic Prices ############
################################################

#######
### Synthetic Prices
#######

syn=fread("[PATH]/estimated prices.csv") # Continuous measure of location-neutral price constructed as described in Paper/Appendix 
syn$sr_property_id=as.character(syn$sr_property_id)

######
### Merge in Syn Prices & Create FE
######

syn=syn[!duplicated(syn),] 

syn=syn[syn$synprice.lvl.allyears>0,] 
syn=syn[syn$synprice.lvl.oneyear>0,]

log.all.high=quantile(syn$synprice.log.allyears,.99)
log.all.low=quantile(syn$synprice.log.allyears,.01)
lvl.all.high=quantile(syn$synprice.lvl.allyears,.99)
lvl.all.low=quantile(syn$synprice.lvl.allyears,.01)
log.one.high=quantile(syn$synprice.log.oneyear,.99)
log.one.low=quantile(syn$synprice.log.oneyear,.01)
lvl.one.high=quantile(syn$synprice.lvl.oneyear,.99)
lvl.one.low=quantile(syn$synprice.lvl.oneyear,.01)

syn=syn[syn$synprice.log.allyears <= log.all.high,]
syn=syn[syn$synprice.log.allyears >= log.all.low,]
syn=syn[syn$synprice.lvl.allyears <= lvl.all.high,]
syn=syn[syn$synprice.lvl.allyears >= lvl.all.low,]
syn=syn[syn$synprice.log.oneyear <= log.one.high,]
syn=syn[syn$synprice.log.oneyear >= log.one.low,]
syn=syn[syn$synprice.lvl.oneyear <= lvl.one.high,]
syn=syn[syn$synprice.lvl.oneyear >= lvl.one.low,]

syngaps=merge(gaps2,syn,by.x=c("sr_property_id","year.int"),by.y=c("sr_property_id","year"))

FEcreate=function(df,num,col) {
  min=min(df[[col]])
  max=max(df[[col]])
  width=(max-min)/num
  zz=cut(df[[col]],seq(min,max,width),include.lowest = T)
  zz=as.factor(zz)
  zz
}

syngaps$FE200.log.all=FEcreate(syngaps,200,"synprice.log.allyears")
syngaps$FE500.log.all=FEcreate(syngaps,500,"synprice.log.allyears")

syngaps$FE200.log.one=FEcreate(syngaps,200,"synprice.log.oneyear")
syngaps$FE500.log.one=FEcreate(syngaps,500,"synprice.log.oneyear")

syngaps$FE200.lvl.all=FEcreate(syngaps,200,"synprice.lvl.allyears")
syngaps$FE500.lvl.all=FEcreate(syngaps,500,"synprice.lvl.allyears")

syngaps$FE200.lvl.one=FEcreate(syngaps,200,"synprice.lvl.oneyear")
syngaps$FE500.lvl.one=FEcreate(syngaps,500,"synprice.lvl.oneyear")


######
### Establish attribute bins 
######

syngaps=merge(syngaps,hedatt,by.x=c("sr_property_id","year.int"),by.y=c("sa_property_id","year"))

syngaps=syngaps[setdiff(1:nrow(syngaps),intersect(which(syngaps$sa_nbr_bedrms==0),which(syngaps$sa_nbr_bath==0))),] 
syngaps=syngaps[which(syngaps$synthrooms>syngaps$sa_nbr_stories),] 
syngaps$sa_yr_blt[which(syngaps$sa_yr_blt==3006)]=2006

syngaps=syngaps[!is.na(syngaps$sa_yr_blt)]
syngaps=syngaps[syngaps$sa_yr_blt!=0,] 

syngaps[, "builtyearbin" := cut(sa_yr_blt,c(0,seq(1900,2020,10)-1),include.lowest = T,labels = LETTERS[1:(length(c(0,seq(1900,2020,10)-1))-1)])] 
syngaps[, "sqftbin" := cut(sa_sqft,c(seq(0,6000,500),seq(6500,10000,1000),seq(10000,50000,5000)),include.lowest = T,
                           labels=letters[1:(length(c(seq(0,6000,500),seq(6500,10000,1000),seq(10000,50000,5000)))-1)])]
syngaps[, "bathbin" := cut(bathadj,c(0,1,2,3,4,5,7,10,15,20))]

syngaps$attbin=paste(syngaps$sqftbin,syngaps$bathbin,syngaps$builtyearbin,
                     syngaps$sa_patio_porch_code,syngaps$sa_pool_code,syngaps$sa_fireplace_code,sep=".")  


######
### Merge in Block Groups
######

syngaps=merge(syngaps,bgs,by.x="sr_property_id",by.y="sa_property_id",all.x=T,all.y=F)
syngaps$BG.ID[which(nchar(syngaps$BG.ID)==0)]=NA  
syngaps$BG.ID=str_pad(syngaps$BG.ID,width = 12,pad = "0",side = "left")

syngaps$govbg=paste(syngaps$group,syngaps$BG.ID,sep=".")
syngaps$govbg=factor(syngaps$govbg)

syngaps$attbin=as.factor(syngaps$attbin)

######
### Run Regressions
######

baseline1=felm(ag2 ~  black.factorvariable | group:year | 0 | group,data=syngaps)
baseline2=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group,data=syngaps)

bin1=felm(ag2 ~ black.factorvariable | group:year + attbin | 0 | group,data=syngaps)
bin2=felm(ag2 ~ blackhisp.factorvariable | group:year + attbin | 0 | group, data=syngaps)

binjur1=felm(ag2 ~ black.factorvariable | group:year + group:attbin | 0 | group,data=syngaps)
binjur2=felm(ag2 ~ blackhisp.factorvariable | group:year + group:attbin | 0 | group, data=syngaps)

bintract1=felm(ag2 ~ black.factorvariable | gtract:year + gtract:attbin | 0 | group,data=syngaps)
bintract2=felm(ag2 ~ blackhisp.factorvariable | gtract:year + gtract:attbin | 0 | group, data=syngaps)

binBG1=felm(ag2 ~ black.factorvariable | govbg:year + govbg:attbin | 0 | group,data=syngaps)
binBG2=felm(ag2 ~ blackhisp.factorvariable | govbg:year + govbg:attbin | 0 | group, data=syngaps)

allbin=list(baseline1,baseline2,bin1,bin2,binjur1,binjur2,bintract1,bintract2,binBG1,binBG2)
save(allbin,file="[OUTPUT PATH]")
gc()

## Repeat process for SYN200FE

syn200FE1=felm(ag2 ~  black.factorvariable | group:year + FE200.log.one | 0 | group,data=syngaps)
syn200FE2=felm(ag2 ~ blackhisp.factorvariable | group:year + FE200.log.one | 0 | group,data=syngaps)

synjur200FE1=felm(ag2 ~  black.factorvariable | group:year + group:FE200.log.one | 0 | group,data=syngaps)
synjur200FE2=felm(ag2 ~ blackhisp.factorvariable | group:year + group:FE200.log.one | 0 | group,data=syngaps)

syntract200FE1=felm(ag2 ~  black.factorvariable | gtract:year + gtract:FE200.log.one | 0 | group,data=syngaps)
syntract200FE2=felm(ag2 ~ blackhisp.factorvariable | gtract:year + gtract:FE200.log.one | 0 | group,data=syngaps)

synBG200FE1=felm(ag2 ~  black.factorvariable | govbg:year + govbg:FE200.log.one | 0 | group,data=syngaps)
synBG200FE2=felm(ag2 ~ blackhisp.factorvariable | govbg:year + govbg:FE200.log.one | 0 | group,data=syngaps)

allsyn200=list(syn200FE1,syn200FE2,synjur200FE1,synjur200FE2,syntract200FE1,syntract200FE2,synBG200FE1,synBG200FE2)
save(allsyn200,file="[OUTPUT PATH]")
gc()

## ...and for SYN500FE
syn500FE1=felm(ag2 ~  black.factorvariable | group:year + FE500.log.one | 0 | group,data=syngaps)
syn500FE2=felm(ag2 ~ blackhisp.factorvariable | group:year + FE500.log.one | 0 | group,data=syngaps)

synjur500FE1=felm(ag2 ~  black.factorvariable | group:year + group:FE500.log.one | 0 | group,data=syngaps)
synjur500FE2=felm(ag2 ~ blackhisp.factorvariable | group:year + group:FE500.log.one | 0 | group,data=syngaps)

syntract500FE1=felm(ag2 ~  black.factorvariable | gtract:year + gtract:FE500.log.one | 0 | group,data=syngaps)
syntract500FE2=felm(ag2 ~ blackhisp.factorvariable | gtract:year + gtract:FE500.log.one | 0 | group,data=syngaps)

synBG500FE1=felm(ag2 ~  black.factorvariable | govbg:year + govbg:FE500.log.one | 0 | group,data=syngaps)
synBG500FE2=felm(ag2 ~ blackhisp.factorvariable | govbg:year + govbg:FE500.log.one | 0 | group,data=syngaps)

allsyn500=list(syn500FE1,syn500FE2,synjur500FE1,synjur500FE2,syntract500FE1,syntract500FE2,synBG500FE1,synBG500FE2)
save(allsyn500,file="[OUTPUT PATH]")
gc()

## Additions 
baseline.tract1=felm(ag2 ~  black.factorvariable | gtract:year | 0 | group,data=syngaps)
baseline.tract2=felm(ag2 ~ blackhisp.factorvariable | gtract:year | 0 | group,data=syngaps)

baseline.bg1=felm(ag2 ~  black.factorvariable | govbg:year | 0 | group,data=syngaps)
baseline.bg2=felm(ag2 ~ blackhisp.factorvariable | govbg:year | 0 | group,data=syngaps)

price.tract1=felm(ag2 ~  black.factorvariable + logmkt | gtract:year | 0 | group,data=syngaps)
price.tract2=felm(ag2 ~ blackhisp.factorvariable + logmkt | gtract:year | 0 | group,data=syngaps)

price.bg1=felm(ag2 ~  black.factorvariable + logmkt | govbg:year | 0 | group,data=syngaps)
price.bg2=felm(ag2 ~ blackhisp.factorvariable + logmkt | govbg:year | 0 | group,data=syngaps)

additions=list(baseline.tract1,baseline.tract2,baseline.bg1,baseline.bg2)
additions.clust=getNclust(additions)
stargazer(additions,type="text", 
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Specification?","baseline-tract","baseline-tract","baseline-bg","baseline-bg",
                           "tract-price","tract-price","bg-price","bg-price"),c("No. Clusters",additions.clust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

#########
### Tables 
#########

specorderbin=c("baseline","baseline",rep(c("attbin","attbin-Jur","attbin-tract","attbin-BG"),each=2))
specorder200=rep(c("200Q","200Q-Jur","200Q-tract","200Q-BG"),each=2)
specorder500=rep(c("500Q","500Q-Jur","500Q-tract","500Q-BG"),each=2)


allbin.clust=getNclust(allbin)
stargazer(allbin,type="latex", 
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Price FE",specorderbin),c("No. Clusters",allbin.clust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)


allsyn200.clust=getNclust(allsyn200)
stargazer(allsyn200,type="latex", 
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Price FE",specorder200),c("No. Clusters",allsyn200.clust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

allsyn500.clust=getNclust(allsyn500)
stargazer(allsyn500,type="latex", 
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,         
          add.lines=list(c("Price FE",specorder500),c("No. Clusters",allsyn500.clust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "log(Assessment Ratio)",
          notes=NULL,
          out = NULL #"[OUTPUT PATH]"
)

################################################
### Test of Racial Diff in Txns ################
################################################


zillow=fread("[PATH]/zillow, simulated values.csv")
zillow$sr_property_id=as.character(zillow$sr_property_id)
zillow=zillow[!is.na(zillow$zillow_val),] 

transgaps=merge(gaps2,zillow,by=c("sr_property_id","sr_date_transfer")) 
transgaps=transgaps[transgaps$salegap>0,]
transgaps$growth=log(transgaps$sr_val_transfer)-log(transgaps$zillow_val)

### Regressions ###

trans1=felm(growth ~ black.factorvariable | govbg:year | 0 | group,data=transgaps)
trans2=felm(growth ~ blackhisp.factorvariable | govbg:year | 0 | group,data=transgaps)

trans.main=list(trans1,trans2)
trans.clust.main=getNclust(trans.main)

stargazer(trans.main,type="text",
          title="Transaction Price Difference",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          order=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          covariate.labels = c("Black Seller","Black or Hispanic Seller","Other Non-White Seller"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jurisd-BG-Yr",length(trans.main))),c("No. Clusters",trans.clust.main)),
          digits=3,digits.extra = 1,float=F,dep.var.caption = "Unexpected Component of Transaction Price",
          out = NULL #"[OUTPUT PATH]"
)


## Produce Appendix Table comparing two samples for robustness
intranstest=paste(transgaps$sr_property_id,transgaps$sr_date_transfer,sep=".")
trans.robust=gaps2
trans.robust$uniqueid=paste(trans.robust$sr_property_id,trans.robust$sr_date_transfer,sep=".")
trans.robust$intranstest=ifelse(trans.robust$uniqueid %in% intranstest,1,0)

rob1=felm(ag2 ~ black.factorvariable | group:year | 0 | group,data=trans.robust[trans.robust$intranstest==1,])
rob2=felm(ag2 ~ black.factorvariable | group:year | 0 | group,data=trans.robust[trans.robust$intranstest==0,])

rob3=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group,data=trans.robust[trans.robust$intranstest==1,])
rob4=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group,data=trans.robust[trans.robust$intranstest==0,])

rob5=felm(black.acs ~ intranstest | group:year | 0 | group,data=trans.robust)
rob6=felm(blackhisp.acs ~ intranstest | group:year | 0 | group,data=trans.robust)

rob=list(rob1,rob2,rob3,rob4,rob5,rob6)
rob.clust=getNclust(rob)
stargazer(rob,type="text",
          title="Transaction Price Difference",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","inamirtest"),
          order=c("black.factorvariableB","blackhisp.factorvariableBH","inamirtest"),
          covariate.labels = c("Black Seller","Black or Hispanic Seller","Used in Txn Price Test"),
          column.labels = c("Assessment Gap","Tract Demographics"), 
          column.separate = c(4,2),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Sample","Test","Not Test","Test","Not Test","Full","Full"),
                         c("Fixed Effects",rep("Jurisd-Yr",length(rob))),c("No. Clusters",rob.clust)),
          digits=3,digits.extra = 1,float=F,dep.var.caption = "",
          out = NULL #"[OUTPUT PATH]"
)

################################################
### Produce Quantile Estimates by Income #######
################################################


incgaps=allgaps[,c("sr_property_id","fips_10","group","year","geoid_10","state","tract","ag2",
                   "hhinc.med.total.acs","sr_date_transfer","medianvalue","sr_val_transfer",
                 "black.factorvariable","blackhisp.factorvariable","allrace.factorvariable","other.factorvariable", "blackflag",  
                 "black.acs","nonwhite.acs","blackhisp.acs","white.nh","pop.total.acs"
)]

incgaps$hhinc.med.total.acs=as.numeric(incgaps$hhinc.med.total.acs)
incgaps=incgaps[incgaps$year>=2005,]         
incgaps=incgaps[complete.cases(incgaps),] 

incgaps$gtract=paste(incgaps$group,incgaps$tract,sep=".") 
incgaps$gtract=as.factor(incgaps$gtract)
incgaps$year=as.factor(incgaps$year)  
incgaps$group=as.factor(incgaps$group)

incgaps$black.factorvariable=as.factor(incgaps$black.factorvariable)
incgaps$blackhisp.factorvariable=as.factor(incgaps$blackhisp.factorvariable)
incgaps$allrace.factorvariable=as.factor(incgaps$allrace.factorvariable)
incgaps$other.factorvariable=as.factor(incgaps$other.factorvariable)

incgaps$black.factorvariable=relevel(incgaps$black.factorvariable,ref="W")
incgaps$blackhisp.factorvariable=relevel(incgaps$blackhisp.factorvariable,ref="W")
incgaps$allrace.factorvariable=relevel(incgaps$allrace.factorvariable,ref="W")
incgaps$other.factorvariable=relevel(incgaps$other.factorvariable,ref="W")


uniqueincometracts=unique(incgaps$hhinc.med.total.acs)


incgaps[,"tractincomebin" := cut(hhinc.med.total.acs,quantile(uniqueincometracts,seq(0,1,1/20)),
                               include.lowest = T,labels = LETTERS[1:20])]

incgaps$tractincvar=incgaps$tractincomebin:incgaps$black.factorvariable
incgaps$tractincvar=relevel(incgaps$tractincvar,ref="A:W")
tractinc=felm(ag2 ~ tractincvar | group:year | 0 | group,data=incgaps)
tractinc.table=samplesplit(tractinc,"tractincvar",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(tractinc.table,names=c("Low",2:19,"High"),clr="blue",
               main="Assessment Gap by Tract Income Vigintile",ylim=c(0,.35),num=20,
               ylab="Assessment Gap, Black Homeowners")
dev.off()


incgaps$tractincvarBH=incgaps$tractincomebin:incgaps$blackhisp.factorvariable
incgaps$tractincvarBH=relevel(incgaps$tractincvarBH,ref="A:W")
tractincBH=felm(ag2 ~ tractincvarBH | group:year | 0 | group,data=incgaps)
tractinc.table.BH=samplesplit.bh(tractincBH,"tractincvarBH",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(tractinc.table.BH,names=c("Low",2:19,"High"),clr="blue",
               main="Assessment Gap by Tract Income Vigintile",ylim=c(0,.35),num=20,
               ylab="Assessment Gap, Black or Hispanic Homeowners")
dev.off()


####################
### Within neighborhood income split
####################

income=fread("[PATH]/HMDA income linked to id & date.csv") # Individual income from HMDA records
income$sr_property_id=as.character(income$sr_property_id)

incgaps2=merge(incgaps,income[,c("sr_property_id","sr_date_transfer","income.mean")],by=c("sr_property_id","sr_date_transfer"),all.x=T,all.y=F)
incgaps2=incgaps2[!is.na(incgaps2$income.mean),]
incgaps2=incgaps2[incgaps2$income.mean<500,]

incgaps2[,"personalincomebin" := cut(income.mean,quantile(income.mean,seq(0,1,1/20)),include.lowest = T,labels = LETTERS[1:20])]
incgaps2$personalincvar=incgaps2$personalincomebin:incgaps2$black.factorvariable
incgaps2$personalincvar=relevel(incgaps2$personalincvar,ref="A:W")

incgaps2=merge(incgaps2,bgs,by.x="sr_property_id",by.y="sa_property_id",all.x=T,all.y=F)
incgaps2$BG.ID[which(nchar(incgaps2$BG.ID)==0)]=NA  
incgaps2$BG.ID=str_pad(incgaps2$BG.ID,width = 12,pad = "0",side = "left")

incgaps2$govbg=paste(incgaps2$group,incgaps2$BG.ID,sep=".")
incgaps2$govbg=factor(incgaps2$govbg)

personalinc.bg=felm(ag2 ~ personalincvar | govbg:year | 0 | group, data=incgaps2)
personalinc.bg.table=samplesplit(personalinc.bg,"personalincvar",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(personalinc.bg.table,names=c("Low",2:19,"High"),clr="blue",
               main="Within Neighborhood by Homeowner Income Vigintile",ylim=c(0,.1),num=20,
               ylab="Assessment Gap, Black Homeowners")
dev.off()


incgaps2$personalincvarBH=incgaps2$personalincomebin:incgaps2$blackhisp.factorvariable
incgaps2$personalincvarBH=relevel(incgaps2$personalincvarBH,ref="A:W")
personalincBH.bg=felm(ag2 ~ personalincvarBH | govbg:year | 0 | group, data=incgaps2)
personalinc.table.BH.bg=samplesplit.bh(personalincBH.bg,"personalincvarBH",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(personalinc.table.BH.bg,names=c("Low",2:19,"High"),clr="blue",
               main="Within Neighborhood by Homeowner Income Vigintile",ylim=c(0,.1),num=20,
               ylab="Assessment Gap, Black or Hispanic Homeowners")
dev.off()


########################
### Split on median tract-level home value
########################

incgaps[,"tracthomebin" := cut(medianvalue,quantile(medianvalue,seq(0,1,1/20)),include.lowest = T,labels = LETTERS[1:20])]
incgaps$tracthomevalue=incgaps$tracthomebin:incgaps$black.factorvariable
incgaps$tracthomevalue=relevel(incgaps$tracthomevalue,ref="A:W")
tractval=felm(ag2 ~ tracthomevalue | group:year | 0 | group,data=incgaps)
tractval.table=samplesplit(tractval,"tracthomevalue",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(tractval.table,names=c("Low",2:19,"High"),clr="blue",
               main="Assessment Gap by Tract Median Home Value Vigintile",ylim=c(0,.25),num=20,
               ylab="Assessment Gap, Black Homeowners")
dev.off()


incgaps$tracthomevalueBH=incgaps$tracthomebin:incgaps$blackhisp.factorvariable
incgaps$tracthomevalueBH=relevel(incgaps$tracthomevalueBH,ref="A:W")
tractvalBH=felm(ag2 ~ tracthomevalueBH | group:year | 0 | group,data=incgaps)
tractvalBH.table=samplesplit.bh(tractvalBH,"tracthomevalueBH",20)

pdf("[OUTPUT PATH]",paper="a4r",width=11,height = 8)
barplot.arrows(tractvalBH.table,names=c("Low",2:19,"High"),clr="blue",
               main="Assessment Gap by Tract Median Home Value Vigintile",ylim=c(0,.25),num=20,
               ylab="Assessment Gap, Black or Hispanic Homeowners")
dev.off()


rm(incgaps,incgaps2,tractincBH,tractinc,personalinc,personalincBH)
gc()

#######################
### Double Portfolio Split
#######################

incgaps[, "obsvaluebin" := cut(sr_val_transfer,quantile(sr_val_transfer,seq(0,1,1/5)),labels=paste("val",LETTERS[1:5],sep=""),include.lowest=T),by="year"] 
incgaps[, "neighvaluebin" := cut(medianvalue,quantile(medianvalue,seq(0,1,1/5)),labels=paste("val",LETTERS[1:5],sep=""),include.lowest=T),by="year"]
incgaps$blacksharebin=ifelse(incgaps$black.acs<.01,"lowblack",
                           ifelse(incgaps$black.acs<.1,"N2",
                                  ifelse(incgaps$black.acs<.25,"N3",
                                         ifelse(incgaps$black.acs<.8,"N4","highblack"
                                         ))))
incgaps[, "blacksharequantile" := cut(black.acs,quantile(black.acs,seq(0,1,1/5)),
                                      labels=c("lowblack","N2","N3","N4","highblack"),
                                      include.lowest=T),by="year"]

incgaps$blackhispsharebin=ifelse(incgaps$blackhisp.acs<.01,"lowblack",
                             ifelse(incgaps$blackhisp.acs<.1,"N2",
                                    ifelse(incgaps$blackhisp.acs<.25,"N3",
                                           ifelse(incgaps$blackhisp.acs<.8,"N4","highblack"
                                           ))))

incgaps$fivebyfive=paste(incgaps$blacksharebin,incgaps$neighvaluebin,sep=".")
incgaps$fivebyfive=as.factor(incgaps$fivebyfive)
incgaps$fivebyfive=relevel(incgaps$fivebyfive,ref="N3.valA")
incgaps$fivebyfive.BH=paste(incgaps$blackhispsharebin,incgaps$neighvaluebin,sep=".")
incgaps$fivebyfive.BH=as.factor(incgaps$fivebyfive.BH)
incgaps$fivebyfive.BH=relevel(incgaps$fivebyfive.BH,ref="N3.valA")


ap2=felm(ag2 ~ fivebyfive | group:year | 0 | group, data=incgaps)

table=data.table("low"=ap2$coefficients[grep("lowblack",rownames(ap2$coefficients),value = F)],
                 "N2"=ap2$coefficients[grep("N2",rownames(ap2$coefficients),value = F)],
                 "N3"=c(0,ap2$coefficients[grep("N3",rownames(ap2$coefficients),value = F)]),
                 "N4"=ap2$coefficients[grep("N4",rownames(ap2$coefficients),value = F)],
                 "high"=ap2$coefficients[grep("highblack",rownames(ap2$coefficients),value = F)]
)

table=table-min(table)
fwrite(table,"[OUTPUT PATH]/double portfolio, B.csv")


###########################################
### Cook County Analysis ##################
###########################################

## CCAO Files
app.indrace=fread("[PATH]/cook county, appeals subset clean, V3.csv",integer64 = "character")
fullcook=fread("[PATH]/cook county, full panel clean, V3.csv",integer64 = "character")

########
### File for win probability
########

first=app.indrace[,c("win","black.factorvariable","blackhisp.factorvariable","other.factorvariable","tract","BG.ID","taxyear","ass_win","bor_win")]
first=first[first$taxyear>=2005,] 
first=first[complete.cases(first),]
first$tract=as.factor(first$tract)
first$taxyear=as.factor(first$taxyear)
first$BG.ID=as.factor(first$BG.ID)

first$black.factorvariable=as.factor(first$black.factorvariable)
first$black.factorvariable=relevel(first$black.factorvariable,ref="W")
first$blackhisp.factorvariable=as.factor(first$blackhisp.factorvariable)
first$blackhisp.factorvariable=relevel(first$blackhisp.factorvariable,ref="W")

first=first[first$black.factorvariable!="",] 

#########
### File for propensity to appeal
#########

full=fullcook[,c("appeal","black.factorvariable","blackhisp.factorvariable","other.factorvariable","tract","BG.ID","taxyear")]

full=full[full$taxyear>=2005,] 
full=full[full$black.factorvariable!="",]

full=full[complete.cases(full),]
full$tract=as.factor(full$tract)
full$BG.ID=as.factor(full$BG.ID)
full$taxyear=as.factor(full$taxyear)

full$black.factorvariable=as.factor(full$black.factorvariable)
full$black.factorvariable=relevel(full$black.factorvariable,ref="W")
full$blackhisp.factorvariable=as.factor(full$blackhisp.factorvariable)
full$blackhisp.factorvariable=relevel(full$blackhisp.factorvariable,ref="W")

########
### File for Success
########

app.indrace$percentreduced=app.indrace$win_reduction/app.indrace$proppose_av
app.indrace$percentreduced.ccao=app.indrace$ass_reduction/app.indrace$proppose_av
app.indrace$percentreduced.bor=app.indrace$bor_reduction/app.indrace$proppose_av

success=app.indrace[app.indrace$win==1,c("win_reduction","black.factorvariable","blackhisp.factorvariable","other.factorvariable","tract",
                                         "taxyear","BG.ID","percentreduced","proppose_av","percentreduced.ccao","percentreduced.bor")]

success=success[success$taxyear>=2005,] 
success=success[success$black.factorvariable!="",]
success$black.factorvariable=as.factor(success$black.factorvariable)
success$black.factorvariable=relevel(success$black.factorvariable,ref="W")
success$blackhisp.factorvariable=as.factor(success$blackhisp.factorvariable)
success$blackhisp.factorvariable=relevel(success$blackhisp.factorvariable,ref="W")

success$tract=factor(success$tract)
success$taxyear=factor(success$taxyear)
success$logreduction=log(success$win_reduction)
success$BG.ID=as.factor(success$BG.ID)

########
### Regressions
########

propen.bg=felm(appeal ~ black.factorvariable | BG.ID:taxyear | 0 | BG.ID,data=full)
propen.bg.bh=felm(appeal ~ blackhisp.factorvariable | BG.ID:taxyear | 0 | BG.ID,data=full)

win.bg=felm(win ~ black.factorvariable | BG.ID:taxyear | 0 | BG.ID, data=first)
win.bg.bh=felm(win ~ blackhisp.factorvariable | BG.ID:taxyear | 0 | BG.ID, data=first)

red.bg=felm(percentreduced ~ black.factorvariable | BG.ID:taxyear | 0 | BG.ID,data=success) 
red.bg.bh=felm(percentreduced ~ blackhisp.factorvariable | BG.ID:taxyear | 0 | BG.ID,data=success) 

########
### Total Effect
########

fullcook1=fullcook
fullcook1=fullcook1[fullcook1$win_reduction>=0 | is.na(fullcook1$win_reduction),] 
fullcook1$newassessments=fullcook1$proppose_av-fullcook1$win_reduction

fullcook1$newassessments[which(is.na(fullcook1$newassessments))]=fullcook1$sa_val_assd[which(is.na(fullcook1$newassessments))]    
fullcook1$proppose_av[which(is.na(fullcook1$proppose_av))]=fullcook1$sa_val_assd[which(is.na(fullcook1$proppose_av))]             
fullcook1=fullcook1[!which(fullcook1$win==0 & fullcook1$proppose_av!=fullcook1$sa_val_assd),] # Remove the 1981 where appealsdata_post != ATTOM value.

fullcook1$logreduction=log(fullcook1$proppose_av/fullcook1$newassessments)
cutpoint_red=quantile(fullcook1$logreduction[which(fullcook1$win==1)],.995)
fullcook1=fullcook1[fullcook1$logreduction < cutpoint_red,] # Trim a small number of very large reductions that look like data errors.

fullcook1$yearint=fullcook1$taxyear
fullcook1$BG.ID=as.factor(fullcook1$BG.ID)
fullcook1$taxyear=as.factor(fullcook1$taxyear)

fullcook1=fullcook1[fullcook1$black.factorvariable %in% c("B","W","O"),]  
fullcook1=fullcook1[fullcook1$blackhisp.factorvariable %in% c("BH","W","O"),]  
fullcook1$black.factorvariable=as.factor(fullcook1$black.factorvariable)
fullcook1$blackhisp.factorvariable=as.factor(fullcook1$blackhisp.factorvariable)
fullcook1$black.factorvariable=relevel(fullcook1$black.factorvariable,ref="W")
fullcook1$blackhisp.factorvariable=relevel(fullcook1$blackhisp.factorvariable,ref="W")
fullcook1=fullcook1[fullcook1$yearint>=2005,] 

totalappealchange1=felm(logreduction ~ black.factorvariable | taxyear | 0 | BG.ID,data=fullcook1)
totalappealchange2=felm(logreduction ~ blackhisp.factorvariable | taxyear | 0 | BG.ID,data=fullcook1)

totalappealchange3=felm(logreduction ~ black.factorvariable | taxyear:BG.ID | 0 | BG.ID,data=fullcook1)
totalappealchange4=felm(logreduction ~ blackhisp.factorvariable | taxyear:BG.ID | 0 | BG.ID,data=fullcook1)


########
### Main Tables 
########

cook=list(cookB.bg,propen.bg,win.bg,red.bg,totalappealchange3)
cookclust=getNclust(cook)

stargazer(cook,type="latex",
          title="Cook County Appeals",
          keep=c("black.factorvariableB"),
          covariate.labels = c("Black Mortgage Holder"),
          apply.coef = times100,
          apply.se = times100,
          keep.stat = c("n","rsq"),
          dep.var.labels = c("Inequality/BG","Appeal   ","Win Appeal    ","  Reduction", "Total Effect"),
          dep.var.labels.include = T,
          add.lines=list(c("Baseline Rate","N/A", "14.6","67.4","12.0","N/A"),c("Fixed Effects",rep("BG-Year",length(cook))),c("No. Clusters",cookclust)),
          digits=3,digits.extra = 1,float=F,
          column.sep.width="8pt",
          dep.var.caption = "Dependent Variable:",
          out = NULL #"[OUTPUT PATH]"
)


cook2=list(cookBH.bg,propen.bg.bh,win.bg.bh,red.bg.bh,totalappealchange4)
cookclust2=getNclust(cook2)

times100=function(x) { x*100}

stargazer(cook2,type="latex",
          title="Cook County Appeals",
          keep=c("blackhisp.factorvariableBH"),
          covariate.labels = c("Black or Hispanic Mortgage Holder"),
          apply.coef = times100,
          apply.se = times100,
          keep.stat = c("n","rsq"),
          dep.var.labels = c("Inequality/BG", "Appeal   ","Win Appeal    ","  Reduction", "Total Effect"),
          dep.var.labels.include = T,
          add.lines=list(c("Baseline Rate", "NA", "14.6","67.4","12.0","N/A"),c("Fixed Effects",rep("BG-Year",length(cook2))),c("No. Clusters",cookclust2)),
          digits=3,digits.extra = 1,float=F,
          column.sep.width="8pt",
          dep.var.caption = "Dependent Variable:",
          out = NULL #"[OUTPUT PATH]"
)


###########################################
### Eff. Rates and Pass Through ###########
###########################################


taxpaid=fread("[PATH]/tax actually paid, clean.csv",integer64 = "character") 
taxpaid$sa_property_id=as.character(taxpaid$sa_property_id)
taxgaps=merge(allgaps,taxpaid,by.x=c("sr_property_id","year"),by.y=c("sa_property_id","year"))
taxgaps$totalexemp=taxgaps$sa_tax_val-taxgaps$sa_tax_val.afterexemp

taxcols=grep("tax_val",colnames(taxpaid),value = T)
for (i in taxcols) {
  addcol=paste("effrate.",i,sep="")
  taxgaps[[addcol]]=log(taxgaps[[i]])-log(taxgaps$sr_val_transfer)
}

taxgaps=taxgaps[taxgaps$hmdasubset==1,]
taxgaps=taxgaps[taxgaps$effrate.sa_tax_val<log(.25),]

# This removes ~100 obs, and is necessary with log specification
taxgaps=taxgaps[taxgaps$sa_tax_val>0,]
taxgaps=taxgaps[which(taxgaps$sr_val_transfer>taxgaps$totalexemp),] 
taxgaps$black.factorvariable=as.factor(taxgaps$black.factorvariable)
taxgaps$black.factorvariable=relevel(taxgaps$black.factorvariable,ref="W")
taxgaps$blackhisp.factorvariable=as.factor(taxgaps$blackhisp.factorvariable)
taxgaps$blackhisp.factorvariable=relevel(taxgaps$blackhisp.factorvariable,ref="W")


taxvar=c("effrate.sa_tax_val","effrate.sa_tax_val.afterexemp")
usevars=c(taxvar,"year","group","ag2",
          "non.wbh","black.acs","blackhisp.acs",
          "black.factorvariable","blackhisp.factorvariable","other.factorvariable",
          "hhinc.med.total.acs","ownerpercent","unemployment","gini","hhfoodstamps")


taxes=taxgaps[,..usevars]
taxes=taxes[complete.cases(taxes),]
taxes$group=as.factor(taxes$group)
taxes$year=as.factor(taxes$year)
taxes=taxes[!is.infinite(taxes$effrate.sa_tax_val.afterexemp),] 


taxbaseB=felm(ag2 ~ black.factorvariable | group:year | 0 | group, data=taxes)
taxbaseBH=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 | group, data=taxes)

t1=felm(effrate.sa_tax_val ~ black.factorvariable | group:year | 0 | group,data=taxes)
t1a=felm(effrate.sa_tax_val.afterexemp ~ black.factorvariable | group:year | 0 | group,data=taxes)

t2=felm(effrate.sa_tax_val ~ blackhisp.factorvariable | group:year | 0 | group,data=taxes)
t2a=felm(effrate.sa_tax_val.afterexemp ~ blackhisp.factorvariable | group:year | 0 | group,data=taxes)

taxregs=list(taxbaseB,t1,t1a,taxbaseBH,t2,t2a)
clusterlength=getNclust(taxregs)

stargazer(taxregs,type="text",
          title="Effective Tax Rate Results",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
              order=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder","Other Nonwhite Mortgage Holder"),
         
          column.labels = rep(c("Assmt. Gap","Before Exemptions","Tax Bill"),2),
         
          apply.coef = times100,
          apply.se = times100,
          keep.stat = c("n","rsq"),
          
          dep.var.labels.include = F,
          add.lines=list(c("Jurisd-Year FE",rep("Y",length(taxregs))),c("No. Clusters",clusterlength)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "Effective Tax Rate - In Sale Year (\\%)",
         
          out = NULL #"[OUTPUT PATH]"
)

######
### Pass through
######


tar1=felm(effrate.sa_tax_val ~ ag2 | group:year | 0 | group,data=taxes)
tar2=felm(effrate.sa_tax_val ~ ag2:blackhisp.factorvariable | group:year | 0 | group,data=taxes)
tar3=felm(effrate.sa_tax_val.afterexemp ~ ag2:blackhisp.factorvariable | group:year | 0 | group,data=taxes)

taxpass=list(tar1,tar2,tar3)
taxpassclust=getNclust(taxpass)

stargazer(taxpass,type="text",
          title="",
          keep=c("ag2","ag2:blackhisp.factorvariableW","ag2:blackhisp.factorvariableBH","ag2:other.factorvariableO"),
          order=c("ag2","ag2:blackhisp.factorvariableW","ag2:blackhisp.factorvariableBH","ag2:other.factorvariableO"),
          covariate.labels = c("All Mortgage Holders"," White Mortgage Holder","Black or Hispanic Mortgage Holder","Other Nonwhite Mortgage Holder"),
          column.labels = c("Before Exemptions","Before Exemptions","Tax Bill"),
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Fixed Effects",rep("Jurisd-Year",3)),c("No. Clusters",taxpassclust)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "Effective Tax Rate - Year of Sale (\\%)",
          out = NULL #"[OUTPUT PATH]"
)

######
### Appendix lag effective tax rate
######

taxvar=c("effrate.sa_tax_val.lag1","effrate.sa_tax_val.afterexemp.lag1")
usevars=c(taxvar,"year","group","ag2",
          "non.wbh","black.acs","blackhisp.acs",
          "black.factorvariable","blackhisp.factorvariable","other.factorvariable",
          "hhinc.med.total.acs","ownerpercent","unemployment","gini","hhfoodstamps")


lagtaxes=taxgaps[,..usevars]  # Note: pulls out of "taxgaps" rather than "allgaps"
lagtaxes=lagtaxes[complete.cases(lagtaxes),]
lagtaxes$group=as.factor(lagtaxes$group)
lagtaxes$year=as.factor(lagtaxes$year)

lagtaxes=lagtaxes[!is.infinite(lagtaxes$effrate.sa_tax_val.lag1)]
lagtaxes=lagtaxes[!is.infinite(lagtaxes$effrate.sa_tax_val.afterexemp.lag1)]


lagt1=felm(effrate.sa_tax_val.lag1 ~ black.factorvariable | group:year | 0 | group,data=lagtaxes)
lagt1a=felm(effrate.sa_tax_val.afterexemp.lag1 ~ black.factorvariable | group:year | 0 | group,data=lagtaxes)

lagt2=felm(effrate.sa_tax_val.lag1 ~ blackhisp.factorvariable | group:year | 0 | group,data=lagtaxes)
lagt2a=felm(effrate.sa_tax_val.afterexemp.lag1 ~ blackhisp.factorvariable | group:year | 0 | group,data=lagtaxes)

lagregs=list(lagt1,lagt1a,lagt2,lagt2a)
clusterlength=getNclust(lagregs)



stargazer(lagregs,type="text",
          title="Effective Tax Rate Results",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          order=c("black.factorvariableB","blackhisp.factorvariableBH","other.factorvariableO"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder","Other Nonwhite Mortgage Holder"),
          column.labels = rep(c("Before Exemptions","Tax Bill"),2),
          apply.coef = times100,
          apply.se = times100,
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Jurisd-Year FE",rep("Y",length(lagregs))),c("No. Clusters",clusterlength)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "Effective Tax Rate - One Year Before Sale (\\%)",
          out = NULL #"[OUTPUT PATH]"
)

###### 
### Appendix - lead effective tax rate

taxvar=c("effrate.sa_tax_val.lead1","effrate.sa_tax_val.afterexemp.lead1")
usevars=c(taxvar,"year","group","ag2",
          "non.wbh","black.acs","blackhisp.acs",
          "black.factorvariable","blackhisp.factorvariable","other.factorvariable",
          "hhinc.med.total.acs","ownerpercent","unemployment","gini","hhfoodstamps")


leadtaxes=taxgaps[,..usevars]  # Note: pulls out of "taxgaps" rather than "allgaps"
leadtaxes=leadtaxes[complete.cases(leadtaxes),]
leadtaxes$group=as.factor(leadtaxes$group)
leadtaxes$year=as.factor(leadtaxes$year)

leadtaxes=leadtaxes[!is.infinite(leadtaxes$effrate.sa_tax_val.lead1)] 
leadtaxes=leadtaxes[!is.infinite(leadtaxes$effrate.sa_tax_val.afterexemp.lead1)] 


leadt1=felm(effrate.sa_tax_val.lead1 ~ black.factorvariable | group:year | 0 | group,data=leadtaxes)
leadt1a=felm(effrate.sa_tax_val.afterexemp.lead1 ~ black.factorvariable | group:year | 0 | group,data=leadtaxes)

leadt2=felm(effrate.sa_tax_val.lead1 ~ blackhisp.factorvariable | group:year | 0 | group,data=leadtaxes)
leadt2a=felm(effrate.sa_tax_val.afterexemp.lead1 ~ blackhisp.factorvariable | group:year | 0 | group,data=leadtaxes)


leadregs=list(leadt1,leadt1a,leadt2,leadt2a)
clusterlength=getNclust(leadregs)


stargazer(leadregs,type="text",
          title="Effective Tax Rate Results",
          keep=c("black.factorvariableB","blackhisp.factorvariableBH"),
          order=c("black.factorvariableB","blackhisp.factorvariableBH"),
          covariate.labels = c("Black Mortgage Holder","Black or Hispanic Mortgage Holder"),
          column.labels = rep(c("Before Exemptions","Tax Bill"),2),
          apply.coef = times100,
          apply.se = times100,
          keep.stat = c("n","rsq"),
          dep.var.labels.include = F,
          add.lines=list(c("Jurisd-Year FE",rep("Y",length(leadregs))),c("No. Clusters",clusterlength)),
          digits=4,digits.extra = 1,float=F,dep.var.caption = "Effective Tax Rate - One Year After Sale (\\%)",
          out = NULL #"[OUTPUT PATH]"
)


###########################################
###########################################
###########################################
###########################################
###########################################
###########################################

